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ABSTRACT 


Flapping-wing propulsion is studied experimentally and numerically. The 
objective of the research is to provide further insight into the aerodynamics of flapping¬ 
wing air vehicles. Experimental work is conducted in the NFS 1.5 m x 1.5 m (5 ft x 5 ft) 
in-draft wind tunnel. A previously constructed long-span flapping-wing model suspended 
by cables is used to approximate the two-dimensional nature of the numerical simulation. 
For this experiment, the model is configured with two wings executing plunge-only 
motion. Thrust is indirectly determined by using a laser rangefinder to measure 
streamwise displacement of the model. Results are compared with previous experimental 
tests. A numerical analysis is conducted using USPOT, a locally developed unsteady 
panel code that models two independently moving airfoils with three degrees of freedom 
and non-linear deforming wakes. Thrust and efficiencies are computed for harmonically 
oscillating airfoils. Direct comparison is made between experimental and numerical 
thrust measurements. 
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I. 


INTRODUCTION 


A. OVERVIEW 

In this paper, flapping-wing propulsion is investigated experimentally and 
numerically. Interest in the field of flapping-wing air vehicles has been spurred by the 
funding of Micro-Air Vehicle (MAV) development. Potential vehicle configurations to 
fulfill the MAV requirement include conventional fixed-wing, flying-wing, rotary-wing 
and flapping-wing designs. Although MAV development has prompted this research, the 
vehicle under consideration is several times larger than MAV specifications. The 
objective of this research is to provide further insight into the aerodynamics of flapping¬ 
wing air vehicles. The ultimate goal is to provide the practical aerodynamic foundation 
necessary for a flapping-wing MAV. This research is limited to the aerodynamics of a 
flapping bi-wing air vehicle, primarily with respect to thrust production. Powerplant, 
structural, and stability and control issues are not considered. 

Experimental work is conducted in the NPS 1.5 m x 1.5 m in-draft wind tunnel. 

A previously constructed long-span flapping-wing model suspended by cables is used to 
approximate the two-dimensional nature of the numerical simulations. A photo of the 
mechanism in the wind tunnel test section is found in Figure 1. For this experiment, the 
model is configured with two wings executing plunge-only motion. Thrust is indirectly 
determined by using a laser rangefinder to measure streamwise displacement of the 
model. Results are compared with previous experimental work. 

A numerical analysis is conducted using USPOT, a locally developed unsteady 
panel code that models two independently moving airfoils with three degrees of freedom 
and nonlinear deforming wakes. Thrust and efficiencies are computed for harmonically 
oscillating airfoils. Direct comparison is made between experimental and numerical 
thrust measurements. 
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Figure 1. Flapping-Wing Mechanism in Wind Tunnel. 


B. BACKGROUND 

In 1992, the Rand Corporation first proposed to the Pentagon the use of Micro-air 
vehicles (MAVs) for aerial observation. The Defense Advanced Research Projects 
Agency (DARPA) became the lead agency for developing this capability. In this role, 
DARPA has defined MAV specifications and funded development. The aircraft are 
defined to be less than 15 cm square and should weigh no more than 0.113 kg. The MAV 
must be capable of "performing a useful military mission" and flying for 20 to 60 minutes 
at a clip of 9 to 18 m/s, “at an affordable cost.” The technology necessary for the 
practical application of these vehicles is rapidly coalescing. 

Due to the small size and low speed of the MAV requirement, low Reynolds 
number flight is encountered. Reynolds numbers less than 100,000 are typical. In the 
current study, the model experiences Reynolds numbers from 0 to 42,600. Since no past 
or current aircraft operate at these Reynolds numbers, simple geometric scaling of larger 
designs is not applicable to a MAV design. 














DARPA requirements do not specify aircraft configuration. Potential candidates 
have taken the form of conventional fixed-wing, flying-wing, rotary-wing and flapping¬ 
wing designs. Low Reynolds Number flight opens a new field of investigation for all 
potential configurations. The current research is aimed at evaluating the aerodynamics of 
a flapping-wing design. Nature demonstrates that flapping wings are efficient at the 
Reynolds numbers of concern, and makes investigation worthwhile. 

C. FLAPPING-WING PROPULSION 

Observations of nature demonstrate the effectiveness of flapping wings in creating 
lift and thrust. The understanding of the aerodynamics of bird and insect flight is limited, 
and achieving nature’s effectiveness in a mechanical device has proven elusive. While 
all potential configurations face hurdles, the flapping-wing configuration faces the 
additional burden of having no previous successes as a practical flying vehicle at any 
scale. Some reasons for the lack of practical success in early aeronautical work with 
flapping wings include a lack of understanding of aerodynamic principles, as well as a 
lack of necessary technologies to implement the ideas. Another primary reason for the 
lack of success at these larger scales is the mechanical limitations of scale, where the 
dynamic forces of large, heavy structures make them unable to support themselves. The 
pursuit of this configuration is prudent in light of the small scale and relatively low 
dynamic forces, the examples given by nature, and previous theoretical work in the field 
that have shown potential efficiencies greater than other configurations. 

An explanation for flapping-wing force production was first observed by Knoller 
[Ref 1] and Betz [Ref 2] in independent studies in 1909 and 1912, respectively. A 
flapping-wing creates an effective angle of attack, cce, as depicted in Figure 2. The 
resulting normal force vector, N, has both lift and thrust components, L and T, 
respectively. Katzmayr [Ref 3] was the first to demonstrate this experimentally, in 1922. 
Refs. 4-6 provide a more detailed history of flapping wing development. 
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Figure 2. Effective Angle of Attack Produced by Flapping-Wing. 

An infinite variety of flapping-wing configurations and dynamics are available for 
analysis. The parameter space is wide, including fixed and flapping airfoils, adjustment 
to the number and relative position of airfoils, different airfoil types, variation of pitch 
and plunge amplitudes and phase angles, changing pivot location and frequency, amongst 
others. The original intent was to explore this vast space with the numerical code and 
apply promising configurations to the experimental model. Delays caused by 
experimental equipment failures, however, subdued this ambitious project. The work has 
therefore concentrated on the evaluation of a single plunge-only configuration. 

The most promising arrangement yielded by previous numerical and experimental 
work by Jones and Platzer [Ref. 4] is the two-airfoil, opposed-plunge configuration. This 
configuration is shown in Figure 3. The lower airfoil is equivalent to the image airfoil in 
ground-effect analysis. This opposed-plunge configuration holds promise in reproducing 
ground-effect behavior when clear of a ground plane. The beneficial effects of ground 
effect are seen put to use by birds, particularly over water. 

When operating near the ground, downwash velocities are decreased. The 
ground presents a boundary condition requiring the normal component of velocity to go 








to zero. This condition alters the streamline pattern around the airfoil. Analytic ground- 
effect analysis places a mirror image about the plane of symmetry. The image vortex 
system is opposite the original system. In order to satisfy the boundary condition, the 
vertical velocities induced by these two image vortex systems will cancel at the plane of 
symmetry. Figure 3 is a graphical representation of this ground-effect system from 
USPOT. In addition to vertical velocity going to zero, the destructive interference of the 
counter-rotating vortex systems diminishes the influence of the trailing vortices on the 
airfoil. The reduced influence, in turn, lessens the induced angle of attack, increasing lift 
and thrust. The configuration under study appropriates the beneficial aspects of ground- 
effect when flying out of ground-effect. In addition, the anti-symmetric motion balances 



Figure 3. Airfoil in Ground Effect 


In order to mimic the 2-dimensional nature of the numerical code as closely as 
possible, a previously constructed long-span flapping-wing mechanism was used in the 
experiments. The mechanism was constructed to support previous research work into 
this field. In addition to providing a promising configuration, the current work provided 
the opportunity to demonstrate repeatability with previous experiments. Repeatability is 
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desired as a means to validate the experimental methods, and to demonstrate that the 
performance of the vehicle will be consistent, something certainly required for a flying 
vehicle. 
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II. NUMERICAL ANALYSIS 


A. PANEL CODE 


Numerical analysis was conducted using the locally developed Unsteady Potential 
flow (USPOT) code. This panel code computes the incompressible flow solution to two 
arbitrarily placed airfoils independently executing three-degree-of-freedom motion. The 
original single-airfoil code (UPOT) was developed by Teng [Ref. 8] and the expanded 
two-airfoil code by Pang [Ref. 9]. The steady portion of the code stems from the work of 
Hess and Smith [Ref. 10], who originally developed the code to analyze steady, inviscid 
flow over an airfoil. The unsteady portion employs the method of Basu and Hancock 
[Ref. 11]. 

In order to obtain the steady solution for the velocity field around the airfoil, the 
velocity potential, <I>, is introduced. The velocity potential is a scalar function in which 
velocity is given by the gradient of O: 


8x. ’ 




V =V<D 

V j 

The velocity potential is defined for irrotational flows only, where: 

VxV =0 


( 1 ) 

( 2 ) 


Flow field velocities are obtained from the velocity potential by differentiating <I> in the 
same direction as the velocity. The airfoil is created by combining several elementary 
potential flows, namely uniform, source, and vortex flow. 

Introducing the definition of the velocity potential to the incompressible 
continuity equation. 


V*F=0 (3) 

yields the Laplace Equation: 

V=4.=0 (o„+4.„.=o) (4) 

Equation 4 is another statement of mass conservation. Note that this equation is a linear 

partial differential equation. As such, elementary flows that are irrotational and 
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incompressible can be combined to form a complicated flow pattern of the same nature 
(irrotational and incompressible). Linearity allows us to superimpose these solutions. 

Superposition of these flows allows the creation of an airfoil of arbitrary shape. 
The code accomplishes this by first dividing the surface of an airfoil into » panels with 
dimensions that are small compared to the airfoil. A source distribution and a vorticity 
distribution are uniformly placed over each panel. The source strength per unit length, qj, 
is constant over a panel, but may vary in strength from panel to panel. This gives n 
unknowns. A single vortex is evenly distributed amongst the panels to provide 
circulation. This vortex strength at each panel, y, is constant. This adds another 
unknown (n+1 unknowns). Boundary conditions are necessary to solve the differential 
equation. Flow tangency is imposed upon the flow, requiring the flow to be tangent to 
the airfoil surface at the center of each panel. This provides n equations. The second 
boundary condition is the Kutta trailing edge condition. The condition is imposed either 
by requiring equal upper and lower surface pressure at the trailing edge, or by requiring 
the tangential velocity on both surfaces to be zero at the trailing edge. This provides one 
equation. There are now n+1 equations and «+l unknowns. The computer can solve this 
system of linearly independent equations simultaneously for the velocity potential. [Ref. 

5] 

In order to account for the unsteady nature of the flow and the associated change 
in circulation, a wake model is used. In this model, the unsteady effects are accounted 
for by application of the theorems of Helmholtz. In the steady model previously 
described, a circulation is induced on the airfoil by the vortex elements. By Helmholtz, 
the total circulation remains constant. Any change in lift entails a change in circulation. 
This circulation change is offset by vorticity shed into the wake. [Ref. 12] 

To ensure the total circulation remains constant in the code, a vorticity panel is 
extended from the trailing edge. This panel’s function is to shed vorticity into the flow, 
counteracting the constantly changing circulation of the pitching/plunging airfoil. The 
wake panel introduces three additional unknowns. These unknowns include the vorticity 
panel’s vortex strength, and the length and orientation of the panel. The strength of the 
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wake panel is equal in magnitude and opposite in direction to the change in circulation 
about the airfoil. The length of the panel varies according to the magnitude of the local 
resultant velocity at the midpoint and the time step. Lastly, the wake panel is orientated 
in the direction of the local velocity at the midpoint of the wake panel. At the end of each 
time step, the vorticity in the wake panel is concentrated into a point vortex, which is then 
released into the flow. The vortex convects downstream, continuing to influence the 
airfoil. [Ref. 5] 

Figure 4 is a representation of the wake model. T is the circulation about the 
airfoil, A is the wake panel length, Yw is the wake panel vortex strength, 0 is the wake 
panel orientation, y and7+I represent panel numbering, qj is the source strength per unit 
length over the designated panel, y is the panel vortex strength, and is the free-stream 
velocity. The subscript k denotes the current time step and ^-1 denotes the previous time 
step. 



Determination of the flow tangency and Kutta condition must be solved 
simultaneously. An iterative procedure must be used to determine the length and 
orientation of the wake panel. The solutions to the potential equations are then used to 
solve for the velocity field about the airfoil. Application of the unsteady incompressible 
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Bernoulli equation gives the pressure distribution about the airfoil. Integrating the 
pressure distribution over the upper and lower surface yields forces in the direction 
normal to the airfoil surface. The net force is broken into components to obtain the 
approximation for lift and thrust. The USPOT code provides results only in a numerical 
output, unlike the single-airfoil code which provides a graphical representation of results. 

Panel codes are subject to limitations. The Laplace equation is a simplified 
version of the Navier-Stokes equation applicable to incompressible flow that neglects 
viscous terms. The equation, therefore, does not take into account separation/stall, or the 
viscous boundary layer effects. This assumption of inviscid flow is most accurate with 
streamlined bodies. The code is only valid for incompressible flow, and is thus only 
useful for low-speed flows, typically flows with a Mach number less than 0.3. The code 
analyzes flow over a 2-D airfoil. No 3-D effects, namely wing-tip vortices, are taken into 
consideration. It is more accurate for predicting lift than drag, since it neglects drag 
induced by wing-tip vortices and skin friction drag. The code compares well with linear 
theory as well as other panel codes. It has been documented extensively with respect to 
linear theory, other numerical methods, and experimental work [Ref. 5]. 

B. NUMERICAL ANALYSIS EQUATIONS OF MOTION 

The generalized airfoil system of USPOT is shown in Figure 5. All lengths are 
scaled and nondimensionalized with respect to the first airfoil chord. Time, t, is 
nondimensionalized with respect to free-stream velocity, Vcc, and chord length, c: 

tV 

r = (5) 

c 

The airfoils are of arbitrary shape and relative position, the second being at coordinates 
X 2 -X 1 , y 2 -y\ from the first. The first airfoil chord length is one, and the second (C 2 ) airfoil 
chord length is a scale factor of the first. The pitching motions of the airfoils are 
independent, and are defined by an initial constant angle of attack, <%, and the change in 
angle of attack, A a: Here, they are shown at angle of attack aj and Oz. The point about 
which each airfoil pivots is relative to its leading edge as a fraction of chord length {xpj, 
Xpi). The frequency of the pitch is entered nondimensionally as the reduced frequency: 



k-^ 


(6) 


where O) is the circular frequency of oscillation, co = Inf. The time-dependent angle of 
attack is given by: 

Qr(r) = OTo + Aacos(^'r) (7) 

The code allows plunging motion in both the x and y coordinate directions. The 
respective plunge amplitudes, hx and hy, are entered as a scale factor of the first airfoil 
chord length. Phase angles, (px and tpy, between the pitch and either plunge motion may 
be entered independently for each airfoil. The motion of the flapping airfoil is sinusoidal 
in the numerical model. The time dependent plunge position of an airfoil is given by: 

= Kx + K cos(^T + 

y{T) = hoy+hyCos{kT + ^^.) 

where hox and hoy are zero for the first airfoil. 
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c. 


LIMITATIONS OF NUMERICAL SIMULATION 


As discussed in the Panel Code section, the numerical simulation is based on the 
Laplace equation. The code is applicable only to incompressible (low-speed) flows, and 
cannot predict viscous effects such as separation/stall, or the viscous boundary layer. No 
3-D effects, such as cross-flow and wing-tip vortices, are accounted for. 

The code is also limited to flapping at reduced frequencies ranging from 
approximately fc=0.1 to 2.0, approximately the limits of reliable code accuracy for this 
configuration. The code’s lower reduced frequency bound is due to singularity strength 
getting very small as k goes to zero. The sum of many very small numbers leads to high 
numerical errors. This is most prominently revealed by the propulsive efficiency, Tjt. 
Both the thrust coefficient and power coefficient become very small with low k (Ct 
Cp k^). Any imprecision in either coefficient will yield questionable results. Since 
thrust efficiency is the ratio of these two very small numbers iTj,= C/Cp), T]t becomes 
inaccurate. The code’s upper reduced frequency bound is dependent upon airfoil 
configuration. As airfoil separation narrows and frequency increases, the code’s wake 
model becomes unstable. The code fails to converge above approximately k=2.0 in this 
close proximity configuration. 

D. NUMERICAL ANALYSIS CONFIGURATION 

For the work conducted, the panel code was set up to model the configuration and 
motion of the experimental model. The experimental model is described in the Flapping- 
Wing Mechanism and Experimental Configuration sections. The configuration was for 
an unswept, constant chord, two-airfoil system (one above the other), with plunge-only 
motion. The associated geometry included mean distance between airfoils, yo=l.4c, 
plunge amplitude, h=0A, and an effective angle of attack, cXe=0°. This opposing plunge 
configuration is depicted in Figure 6. Symmetric NACA 0014 airfoil sections were 
specified, flapping at reduced frequencies ranging from ^0.1 to 2.0. This limitation on 
reduced frequency is discussed in the previous section. 
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The panel code’s input and output are nondimensional. In order to simulate the 
experimental work, it was necessary to nondimensionalize the experimental data, run the 
panel code, then dimensionalize the data for analysis. The experimental configuration 
above provides the input data necessary to run the panel code, most notably the reduced 
frequency, k, and plunge amplitude, hy. The panel code output consists of the lift, drag, 
and moment coefficients, angle of attack, and the plunge position of each airfoil. Thrust 
is obtained from the panel code output by dimensionalizing the drag coefficient, 

T = -D^^pJV^SC^ (10) 

where T is thrust, D is drag, poo is free-stream air density, S is wing area, and Cd is the 
drag coefficient. The wing area must be multiplied by the number of airfoils, so S=2bc, 
where 2 is the number of airfoils, b is the useful lifting span of the experimental 
mechanism (1200 mm), and c is the chord (64 mm). The nondimensional drag 
coefficient Cd is also based on the useful lifting span, so that Co^Cdb. 

The reduced frequency input to the panel code contains both velocity and 
frequency data. Expanding and rearranging the equation for reduced frequency, the 
following equation for velocity is obtained: 
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( 11 ) 


Dimensional velocities are obtained by entering experimental flapping frequencies (2.95, 
5.01, 7.03 Hz), chord length, and reduced frequency. These procedures dimensionalized 
thrust, velocity, and frequency for analysis. 


E. NUMERICAL ANALYSIS RESULTS 

Jones conducted an evaluation of several configurations. He compared the 
numerical code with Garrick’s linear theory for a single flapping airfoil, and conducted a 
numerical analysis for the two-airfoil flapping case. The results are dimensionalized for 
analysis based upon the experimental flapping mechanism’s geometry and dynamics 
described in the previous section. Figure 7 is a plot of thrust coefficient, C, versus 
reduced frequency, k. The graph represents thrust per airfoil. The graph indicates good 
agreement between linear theory and the code for the single airfoil case. The results 
show that for the two-airfoil configuration, the thrust per airfoil is significantly higher 
than for the single airfoil case. This suggests a beneficial thrust aspect associated with 
the opposed-plunge motion. Since the graph represents thrust per airfoil, the two-airfoil 
configuration should produce well over twice the thrust of the single airfoil. As 
discussed in the Flapping-Wing Propulsion section, the opposed-plunge motion replicates 
the motion used in the analysis of airfoils in ground-effect. 
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In Figure 8, the propulsive efficiency, rit, is plotted against k for the same 
configurations. The two-airfoil, opposed-plunge configuration shows higher efficiencies 
than the single flapping airfoil over nearly all of the reduced frequency range. The rapid 
drop in efficiency by the two-airfoil configuration at low k is most likely due to code 
limitations than actual behavior. A discussion of this limitation is contained in the 
Limitations of Numerical Analysis section. 
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Reduced Frequency, k 

Figure 8. Propulsive Efficiency vs Reduced Frequency. After Ref 5. 


The results of USPOT dimensionalized for the experimental model is presented in 
Figure 9. Tabular data is contained in Appendix B. Lines of constant reduced frequency 
and constant effective angle of attack are coincident, and are added for reference. For a 
single flapping airfoil, linear theory predicts thrust will increase roughly proportional to 
f. For a two-airfoil system, we see a similar behavior. 
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Figure 9. Thrust vs Velocity from Panel Code. 
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III. EXPERIMENTAL ANALYSIS 


A. EXPERIMENTAL APPARATUS 


1. Flapping-Wing Mechanism 


The flapping-wing mechanism used for the experiments was designed by Dr. 
Kevin Jones at the Naval Postgraduate School in 1997. Its construction and dynamics are 
described in Ref. 5, from which this description is largely taken. Isometric, profile, and 
top views are shown in Figure 10, Figure 11 and Figure 12. The structure of the fuselage 
is made of aluminum, the front and rear nacelles of balsa-wood, and the wings of wood. 
The device is configured to flap one or two wings. When flapping two wings, the wings 
are configured one above the other. The two flapping wings are configured to execute 
two degree of freedom motion, with both pitch and plunge amplitudes adjustable. The 
phase angle between pitch and plunge is also adjustable. Stationary wings may be 
attached before or after the flapping wings in the slots depicted. 


notch in rear nacelle 



flapping wings 


Figure 10. Isometric View of Upturned Happing Mechanism. From Ref. 5. 
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actuator bearinged flapping flapping motor 
wheels pins wings beams 

Figure 11. Profile View of Flapping Mechanism. From Ref. 5. 


front 



Figure 12. Top View of Flapping Mechanism. After Ref. 5. 


The flapping wings are attached to moving beams, which are driven by eight 
actuator wheels. A pin attached to each of the eight wheels rides in a slot at the end of 
each beam. This attachment changes the circular motion of the wheels into the flapping 
motion of the beams. The beams may be attached to the wheels at a selected wheel 
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radius to control plunge amplitude. Wing pitch is controlled by adjusting the phase 
between the fore and aft wheels. 

The eight wheels are attached to four worm gears. The gears mesh with two 
worms that are driven by an Astro-Flight Cobalt 40 model airplane motor. The motor is 
powered by a variable voltage/current power supply. After burning up one motor 
armature, the new motor’s voltage was limited to 26V for a maximum of one minute. 
This voltage limit restricted maximum flapping frequency to 7 Hz. This corresponds to a 
motor rpm of 21,000 with a 50:1 gear reduction. Current draw was always less than 5A. 
New worms and worm gears are being ordered to obtain higher frequencies. 

Two wing sections were available for the experiment. A symmetric cross-section 
approximating a NACA 0014 and a cambered cross-section approximating a NACA 
23012. Further configuration details can be found in the Experimental Configuration 
section. The wings are manufactured by the Miniature Aircraft Supply model helicopter 
company. The wings are made of balsa wood and laminated hardwood, and covered with 
a thin, smooth plastic sheet. The wings can be bolted to the flapping beams, or attached 
as stationary wings in the slots before and after the flapping wings. The stationary wing 
feature is useful for experiments based upon Schmidt’s wave propeller [Ref. 13]. This 
work suggests that a stationary airfoil in the wake of flapping airfoils may increase thmst 
production. 


2. Wind Tunnel 

Experiments were performed in the NFS low-speed in-draft smoke tunnel 
depicted in Figure 13. The tunnel was modeled after the smoke tunnel of the Naval Air 
Engineering Lab in Philadelphia [Ref. 14]. Air is ingested from inside the building 
through a square 4.5 x 4.5 m inlet, converging through a 9:1 bell shaped contraction to a 
1.5 X 1.5 m test section. Tunnel speed is controlled by a variable pitch fan driven by a 
constant speed electric motor. Motor and fan vibration are isolated from the test section 
by rubber sleeves on each side of the motor/fan assembly. The tunnel velocity range is 0 
to 9.5 m/s. 
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Figure 13: NFS Low-Speed In-draft Wind Tunnel. From Ref. 5. 

The initial experimental effort was delayed by a tunnel fan failure. The effort 
required to fix the fan is described in Appendix A. A portion of the effort consisted of 
resetting the variable pitch fan blades to provide tunnel velocities from 0 to 9.5 m/s. This 
will support future very low-speed work. In addition to the fan repair, old fixtures 
supporting old experiments were all removed, the tunnel floor refinished, all seams 
smoothed and holes sealed. The old pitot-static tube and the associated fluid manometer 
were removed and replaced by a better positioned pitot-static tube with a higher accuracy 
differential pressure transducer. A tunnel calibration was performed following tunnel 
repairs. The test section turbulence intensity level remained below 1.75% with tunnel 
velocities over 1.5 m/s. Tunnel calibration is described in depth in Appendix D. 

3. Laser Doppler Velocimetry (LDV) 

Laser Doppler Velocimetry (LDV) is a precise, accurate, non-intrusive method of 
measuring the velocity of fluid flow. LDV is an optical method that utilizes the Doppler 
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principle, measuring the frequency shift of laser light scattered by a particle in the flow. 

A dual-beam, on-axis, backward-scatter system was used to measure the air velocity and 
turbulence intensity in the NPS low-speed, open circuit smoke tunnel. A general 
description is provided here. Details of LDV position and sampling may be found under 
Experimental Data Collection, a detailed description of the LDV system and concepts in 
Appendix C, and LDV use for tunnel calibration in Appendix D. 

In the dual-beam arrangement utilized, a beam of laser light is split in two, the 
frequency of one of the beams is shifted, and both beams transmitted such that they cross 
each other. Particles moving through the volume in space created by crossing the 
unshifted and shifted beams (an ellipsoid) scatter the laser light. The reflected light is 
collected at the probe (on-axis, backward scatter), and sent via optical receiving fiber to a 
Photo Multiplier Tube (PMT). A photo detector in the PMT separates the light by color, 
and converts it into a voltage signal at the Doppler frequency. The signal produced by 
the photo detector is sent to the IFA 755 signal processor via a frequency downshifter and 
Colorlink. The signal processor determines frequency from the electrical signal and 
sends these signals to the coihputer. The computer converts the Doppler frequency shift 
to velocity of the particles. Flow Information Display (FIND) software was used to 
display and store LDV system outputs. A detailed account of the system equipment, 

LDV concepts, and a system diagram can be found in Appendix C. 

LDV has several advantages over other quantitative flow measurement 
techniques. It is a linear instrument, non-intrusive, requires no calibration, is independent 
of fluid density (suitable for a variety of fluids), and measures only the selected 
component of velocity. LDV disadvantages include high cost and sophistication, a 
requirement for an optical view of the flow, and the dependence on the presence of 
particles to scatter the light. [Ref. 15] 

The initial purpose of setting up the LDV system was to determine the flow field 
dynamics over the flapping wings. Determination of the onset of flow separation was the 
primary goal. Delays caused by failures of associated and unassociated equipment, and a 
delay in the purchase of a longer focal length lens, prevented its use in this fashion. It is 
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intended that future work will use the LDV for this purpose. In the current work, the 
LDV was used for tunnel calibration following tunnel repairs and improvements, and also 
for measuring free-stream velocity at the model during testing. Tunnel calibration 
utilizing the LDV is discussed in Appendix D. 

For experimental testing, the LDV probe was mounted on a 3-axis traverse 
mechanism outside the tunnel at the test section window. A photograph of the probe on 
its 3-axis traverse is contained in Appendix C. The LDV probe is fitted with a 350 mm 
lens, limiting the velocity measurements to less than a quarter of the tunnel width. A 750 
mm focal length lens will be acquired soon, allowing LDV measurements to the center of 
the tunnel. 


4. Miscellaneous Equipment 

A Rosco model 4500 smoke generator was used to provide particles for light 
scattering. A smoke rake was adjusted to direct smoke particles through the probe 
volume, ensuring the greatest scattering of laser light. 

A pitot-static tube located upstream of the model was used for general 
adjustments to tunnel velocity. The pitot-static tube was connected to a MKS Baratron 
type 223B differential pressure transducer. The transducer voltage output is linear with 
pressure. LDV velocity measurements were used for fine tunnel speed adjustments prior 
to experimental runs. Data reduction was based solely on LDV velocities. 

A frequency strobe light was used to get the flapping-wing frequency near that 
desired. The strobe was set to the desired frequency, and the wing frequency adjusted to 
match. The accuracy of the strobe was checked using a Monarch Instrument Tach FV 
Digital Optical Tachometer. Strobe error varied with frequency to a maximum of 1.5 
percent in the frequencies of interest. Actual frequency used for data reduction was 
detennined after the experiment using data from a micro-switch inside the model. 

A Matsushita Electric Works Model ANL1651 AC laser rangefinder was used to 
determine model displacement when flapping. The rangefinder was positioned behind 
the model at its center point distance of 13 cm. Its measurable range is 8 - 18 cm, with 
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accuracy decreasing with increased distance from the center point. The output of the 
rangefinder is 1 V/cm displacement. 

A Link Instmments DSO-2102 digital storage oscilloscope was used to measure 
voltage. The DSO-2102 was connected to a desktop PC for display and recording of the 
voltage signals. The oscilloscope simultaneously records two data channels. One 
channel was designated for the laser rangefinder, and this data was subsequently used to 
determine thrust. The second channel was used to collect micro-switch data for eventual 
determination of flapping frequency. The oscilloscope features 32 kilobyte buffers per 
channel. A 1000 Hz sample rate was used, for a run time of approximately 32 seconds 
per thrust record. 

B. EXPERIMENTAL CONFIGURATION 

In the experiments conducted, the mechanism was bi-wing configured with 
symmetric wings approximating NACA 0014 sections. The wings were of equal chord 
length, separated by a mean distance of 1.4 chord lengths. The wings spanned 1270 mm; 
the long-span intended to replicate the two-dimensional nature of the panel code as 
closely as feasible. The flapping mechanism was limited to plunge-only motion. The 
plunge amplitude for both wings was 25.4 mm, or approximately 40 percent of the chord 
length (h=0.4). The experimental configuration is shown in Table 1. 
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Wing Configuration 


Bi-wing flapping, no stationary wings 

Wing Type 


Symmetric, approximating NACA 0014 

Wing Span 

b* 

1270 mm 

Fuselage Width 


70 mm 

Useful Lifting Span 

b 

1200 mm 

Wing Chord 

c 

64 mm 

Aspect Ratio 

AR=b^/S=b/c 

19.84 full span, 18.75 useful span 

Taper Ratio 

A»—c^/C q 

1 

Sweep Angle 

A 

0° 

Mean Wing Separation 


1.4c 

Plunge Amplitude 

h 

h = 0.4 = 25.4 mm 

Pitching Motion 

Aa 

None, (Xgeorn ” ^ 

Flapping Frequency 


3, 5, 7 Hz (Actual 2.95, 5.01, 7.03 Hz) 

Velocity 

Voo 

0 - 9.5 m/s 

Reduced Frequency 

k 

0.12-oo 


Table 1. Experimental Configuration. 


C. EXPERIMENTAL DATA COLLECTION 

The tunnel velocity, thrust, and flapping frequency were needed in order to take 
useful measurements. First, the tunnel velocity was set with the pitot-static tube, then 
fine tuned with the LDV. The pitot tube was located 2.5 meters upstream of the model. 
Immediately prior to taking thrust measurements, LDV free-stream velocity sampling 
was conducted at the model. The LDV was positioned at the test section window at the 
nose of the model. The focal point was 0.35 m into the tunnel (All LDV measurements 
were taken with a 350 mm focal length lens positioned as close as possible to the tunnel 
glass), and 0.25 meters above the model. Runs were conducted every 0.5 m/s from 0 to 5 
m/s, every 1.0 m/s from 5 to 9 m/s, and a final run was made at maximum tunnel speed of 
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9.5 m/s. Each measurement run consisted of 10,000 samples taken over a period of 3 to 
10 seconds, dependent upon smoke density. This LDV free-stream velocity and 
associated velocity fluctuations were recorded for subsequent data reduction. 

The experimental arrangement for measuring thrust is depicted in Figure 14. A 
photo of the model in the tunnel is shown in . The model is suspended with four thin 
cables 112.8 cm in length. The cables allow displacement along the streamwise axis, but 
hinder motion in other axes. In order to determine thrust, the model displacement due to 
flapping was required. Drag induced by tunnel velocity deflected the model rearward. 
Flapping the wings of the model propelled the model forward. The net displacement was 
measured by the laser rangefinder, positioned on a two-axis traverse behind the model, 
reflecting the beam from a notch in the model’s rear nacelle. 



Figure 14. Flapping Mechanism in Wind Tunnel Test Section. From Ref. 5. 

After setting tunnel velocity, the laser rangefinder was positioned at its zero 
voltage center point, 13 cm behind the model. Since the model’s rearward displacement 
due to drag increased with velocity, the rangefinder was reset to its center point at each 
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tunnel velocity. The model’s drag varied slightly dependent on wing position. The 
oscilloscope displayed a clear decrease in drag when the wings were closer. The model 
was flapped at an extremely low frequency, on the order of 0.1 Hz. This displacement, 
essentially representing its no-flap (0 Hz) condition, was recorded. The mean 
displacement of the model, therefore, took into account the variation in drag with wing 
position. The model was then flapped at 3, 5 and 7 Hz. Displacement voltages were 
again recorded. The mean displacement of the no-flap condition (its ‘static’ position) 
was subtracted from flapping displacements prior to calculating thrust at 3, 5 and 7 Hz. 

In order to obtain the desired frequency, the strobe light was set to 3, 5, or 7 Hz, 
and the motor voltage was adjusted until the beams flapping the wings appeared to be 
stationary. This method was used to set frequency during experimental runs. For the 
ensuing data reduction, determination of flapping frequency was based on a voltage 
signal from a micro-switch within the model. The micro-switch was triggered once per 
cycle. The reduced frequencies corresponding to tunnel velocities from 0 - 9.5 m/s and 
flapping frequencies from 2.95 - 7.03 Hz were ^0.13 -«». 

The flapping mechanism was not intended to be a flying vehicle, and weighs 
approximately 4 kg. The maximum model displacement from flapping, therefore, was a 
modest 1.10 cm. The laser rangefinder voltage was recorded concurrently with the 
frequency micro-switch voltage by the digital oscilloscope. The rangefinder and 
frequency micro-switch data were recorded with a sample rate of 1000 Hz for a period of 
approximately 32 seconds. This ensured a good census of thrust data, measuring between 
94 and 225 cycles per run. 

D. EXPERIMENTAL DATA REDUCTION 

After completing data collection, each of the 16 velocity-frequency combinations 
had an associated velocity files from the LDV and a voltage file from the oscilloscope. 
For each of these experimental runs, the LDV FIND software calculated the mean and 
standard deviation of the velocity. LDV velocity data can be found in Appendix G. 

Before converting the oscilloscope voltage to force, a thrust calibration of the 
model was necessary. A thread attached to the model was strung over a pulley, and 
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known weights were hung on the thread. The model was deflected and allowed to return 
to an equilibrium position. The tunnel was closed to prevent air circulation from 
affecting the model’s displacement. A 60-second voltage sample of the model’s 
displacement was recorded by a digital voltmeter. This was conducted three times for 
each weight. The mean and standard deviations of the voltages from the three runs were 
averaged, and a calibration curve plotted. An approximately linear relationship was 
established between the weight and the laser rangefinder displacement voltage. The data 
and calibration curve can be found in Appendix E. 

The voltage files from the oscilloscope consisted of voltages from both the 
rangefinder and frequency micro-switch. A MATLAB program, found in Appendix F, 
was written to convert the data into useful form. The program accomplishes this by first 
converting the data to voltage, then correcting for oscilloscope zero voltage offset. The 
thrust voltage is converted to force by means of the thrust calibration. The program 
identifies each cycle by the trip switch data. Figure 15 is a pictorial representation of the 
program’s cycle counting. The voltage, thrust, and frequency are calculated for each 
cycle. The minimum, maximum, and average deviation of the frequency is calculated 
from the cycle data, as well as mean thrust. Errors associated with measuring thrust are 
also calculated, and are described in the Error Analysis section. Partial cycle data at the 
beginning and end of each thrust record is disregarded. Data is stored to two files, one 
containing all data for a given frequency-velocity combination, and one to a file 
containing a summary of all frequency-velocity combinations. 



Figure 15. Pictorial Representation of Program Cycle Counting. 
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E. 


EXPERIMENTAL RESULTS AND ANALYSIS 


1. Experimental Results 

The results of the opposing plunge experiments are presented in Figure 16. 
Tabular data can be found in Appendix G. The trends imitate the numerical results 
presented in Figure 9. Frequencies are mean frequencies, all cycles varying less than 
0.05 Hz from the mean, with an average frequency error of 0.9%. Error bars are added 
for both thrust and velocity. The sources of error are analyzed in the Error Analysis 
section. 



01 2 3456789 10 

Velocity [m/s] 

Figure 16. Experimental Thmst vs Velocity. 


2. Error Analysis 

There are four primary sources of error associated with the experimental thrust 
measurements. These four include the average deviation error, laser rangefinder error, 
oscilloscope error, and weight/thmst calibration error. The sum of these four errors is 
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presented as error bars in Figure 16. The individual and total thrust errors are presented 
in Appendix G. The velocity measurements taken by the LDV are very accurate. Errors 
associated with LDV measurements are normally very small. The velocity error bars 
displayed in Figure 16 represent one standard deviation added and subtracted from the 
mean LDV velocity. The LDV data is contained in Appendix G. 

The first source of thrust error is thrust average deviation. The mean thrust for 
every complete flapping cycle of the model was calculated. This cycle mean thrust of the 
model varies from cycle to cycle. This variation is due to a variety of factors. The most 
significant contributor is most likely that the tunnel intake is indoors. Although the 
tunnel room is vented to the outside by several banks of windows, the opening and 
closing of doors and other windows in the building caused slight velocity drift which 
affects model displacement and thrust. Average deviation is similar to standard deviation 
in that it provides a measure of the width of the data about the mean. It is a more robust 
estimator, and is given by the equation below [Ref. 16]. 

1 w 

ADev{x^ ...x ^) = —2 k;•” (^2) 

For this application, x represents the mean of each cycle, y is the cycle counter, N is the 
total number of cycles, and x is the mean of all cycles. The average deviation was 
calculated at every frequency and velocity and is considered here to be an estimated error 
bound. It is the predominant source of error for velocities above approximately 4 m/s. 

The second source of thrust error was the laser rangefinder. Rangefinder error 
consists of a fixed error at the center point (13 cm). A linear error proportional to the 
distance from the center point is added to the fixed error at all other points. For the 
ranges encountered here, the rangefinder accuracy is given by the following equation: 

Range Error- O.hnm ± 0.2% * Distance from Center Point (13) 

The thrust error associated with the rangefinder is the most significant source of error for 
velocities below approximately 4 m/s. 

The third source of thrust error was due to the digital oscilloscope. An analysis of 
the 8-bit oscilloscope accuracy was conducted against a calibrated voltage source and 
against three voltmeters (voltmeters having higher voltage measurement accuracy). A 
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zero voltage offset in the oscilloscope was subtracted from all voltage measurements. 
Also, an error due to the maximum resolution of the oscilloscope was added to all voltage 
measurements. This error varies with the selected maximum value of the displayed 
vertical axis. The highest accuracy is obtained when the voltage signal approaches the 
maximum displayable, but does not exceed the display window. The oscilloscope 
accuracy analysis indicated that the actual oscilloscope error was always less than half of 
its resolution. Statistically, this is reasonable in that errors would tend to occur on both 
sides of the actual voltage, offsetting each other. In the interest of a conservative 
analysis, however, the maximum resolution was used as the error. 

The fourth and final primary source of error was that associated with the thrust 
calibration. As discussed in the Experimental Data Reduction section, three 60-sec 
voltage runs were conducted for each weight. The mean voltages were averaged, and a 
near-linear calibration curve plotted. The standard deviations were also calculated. The 
standard deviation at a particular thrust level was applied as an error due to thrust 
calibration. The error was negligible and had the least impact on total error. The 
calibration data and curve can be found in Appendix E. The error caused by the 
calibration can be found in Appendix G. 

3. Comparison with Previous Experimental Work 

In 1998, Jones used the same model in several configurations for experimental 
work. Jones conducted the experiments at 2,4, 6, and 8 Hz, and the results of the work 
are presented in Ref. 5. The plunge-only experimental work is compared with the current 
results in Figure 17. Jones’ work (even frequencies) is presented as solid characters, the 
current work (odd frequencies) with open characters. Jones’ measurements were taken 
on two separate occasions and the data combined, hence the two points for each 
frequency in the static case. 
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Figure 17. Thrust vs Velocity Comparison w/ Jones. 


Linear theory and the numerical code imply that thrust will increase with velocity. 
The experimental data reflects this trend throughout the velocities tested. Linear theory 
and the numerical code also imply that thrust will increase as the frequency squared. 
Individually, the experimental data sets from both the current work and Jones confirmed 
this expectation. When combined, both data sets confirmed this expected trend in the 
static case. The combined data sets, however, did not reflect this at velocity. The current 
work showed higher thrust levels at 3, 5, and 7 Hz than obtained by Jones at 4,6 and 8 
Hz, respectively. One of the primary reasons for conducting the experiments was to 
demonstrate the repeatability of the experiment. Repeatability is necessary in order to 
demonstrate that the output will be consistently obtained, something certainly required 
for a flying vehicle. In this case, the experimental data does not initially confirm 
repeatability. 


33 




An examination of the experimental data indicates that the tunnel velocities 
recorded by the fluid manometer in the earlier work were inaccurate. Replicating the 
measurements is not possible since both the pitot tube and manometer have since been 
replaced. Research into work conducted after the pitot tube and manometer were 
replaced, but before the tunnel velocity range was changed during a fan overhaul, found 
only one set of experiments was conducted. This set of experiments used a newly 
installed LDV system, but measured at a different tunnel position. An assessment was 
made of the maximum velocities recorded in each experiment, and also the position of 
the velocity measurements and the model. This assessment indicated that the recorded 
velocities were roughly 30% too high. To account for this velocity error, a scaling error 
was applied to the manometer data, shifting the non-static, even frequency thrust data to 
lower velocities. This error correction is applied in Figure 18. This correction indicates 
that the thrust does indeed increase with frequency, as anticipated. The correction 
suggests that the experiment is repeatable, and that a steady thrust can be consistently 
produced. 



Corrected Velocity [m/s] 

Figure 18. Thrust vs Corrected Velocity Comparison w/ Jones 
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F. 


COMPARISON OF EXPERIMENTAL AND NUMERICAL RESULTS 


The experimental and numerical results of the plunge-only configured flapping¬ 
wing mechanism at the frequencies analyzed are shown in Figure 19. 
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Figure 19. Thrust vs Velocity for Experimental and Numerical Data. 

The numerical code is unable to predict static thrust. The numerical predictions 
are limited to reduced frequencies from k=0 .1 to 2.0. These restrictions on k are 
discussed in the Limitations of Numerical Simulation section. As expected, the panel 
code predicts thmst consistently higher than that achieved in the experiments. The 
steady-state drag was subtracted from both the numerical and experimental results. The 
computed thrusts, therefore, eliminate the viscous boundary layer drag effects. In 
general, the experimental data follows the numerical data trends. 

Three potential factors have been identified for the discrepancies between 
experimental and numerical results. The numerical model does not account for 3-D 
losses, namely wing-tip vortices. The numerical code also does not account for 
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separation. Figure 9 shows the high effective angles of attack seen by the wings at high 
k. While the NACA 0014 has a relatively high stall angle, 13°, the wings experience an 
effective angle of attack of 38.7° at k=0.1. Although separation is delayed by the 
plunging motion, localized separation is most likely occurring just above the static stall 
angle. Lastly, the model wings are subject to flexing, with the tips experiencing higher 
amplitudes than the roots. The amplitude varies with frequency, and there is an 
associated phase lag. This effect will tend to raise the experimental values artificially 
high, aggravating the difference between the experimental and numerical thrusts. The 
wings may also flex in torsion, but measurements were not attempted. In general, the 
difference between numerical and experimental thrust may be attributed primarily to 
separation, three-dimensional effects and the flex in the airfoils of the experimental 
model. 
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IV. CONCLUSIONS 


Observations of nature demonstrate the effectiveness of flapping wings in creating 
lift and thrust. The relatively recent interest in Micro-Air Vehicle development has 
spurred interest in developing a flapping-wing MAV. Due to the small size and low 
speed of MAVs, all potential MAV configurations encounter low Reynolds number 
flight, a new field of study. The flapping-wing MAV has the additional burden of having 
no previous successes as a useful flying vehicle at any scale. The basic science and 
technologies necessary to employ a mechanical device in the MAV role are rapidly 
maturing. Nature hints that flapping wings are efficient at the Reynolds numbers of 
concern. The potential efficiency of this design makes investigation worthwhile. 

A numerical and experimental investigation of flapping-wing propulsion was 
conducted. Although prompted by MAV development, the vehicle under consideration is 
several times larger than MAV specifications. The size was dictated in the interest of 
gaining an initial understanding of flapping-wing propulsion and low Reynolds number 
flight, verifying a numerical code, eventual smoke visualization, and general ease of 
assembly. The analysis was limited to the aerodynamics of flapping-wing propulsion. 

The numerical analysis was conducted using USPOT, a locally developed panel 
code. The code models two independently moving airfoils with three degrees of 
freedom. Thrust and efficiencies were computed for two harmonically oscillating airfoils 
executing plunge-only motion. The numerical analysis established the general 
performance trends for opposed-plunge flapping-wing thrust. The code confirmed the 
expectations of linear theory, indicating thrust increases roughly as the square of the 
frequency. 

Experimental work was conducted in the NPS 1.5 m x 1.5 m in-draft wind tunnel. 
A previously constructed long-span flapping-wing model suspended by cables was used 
to approximate the two-dimensional nature of the numerical flow. For this experiment, 
the model was arranged in a bi-wing configuration with the wings executing plunge-only 
motion. Thrust was indirectly determined by measuring the streamwise displacement of 
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the model due to flapping. The steady-state drag was subtracted from the computed 
thrust in order to minimize viscous boundary layer drag effects. 

The experimental results were compared to Jones’ previous work with the model 
in the same configuration. The current work was conducted at different frequencies. The 
data sets reflected an increase of thrust with velocity, a trend indicated by both the linear 
theory and the numerical code. Individually, the experimental data sets from both the 
current work and Jones also confirmed the expected trend of increasing thmst with 
increasing velocity. When combined, however, the current work showed higher thrust 
levels at 3, 5, and 7 Hz than obtained by Jones at 2,4, and 6 Hz, respectively. An 
examination of the previous work indicated that the recorded tunnel velocities were too 
high. Shifting the old data to account for this velocity error corrected this discrepancy. 
The corrected data reflects that thrust does increase with frequency, as anticipated. The 
correction suggests that the experiment is repeatable, and that a steady thmst can be 
consistently produced. 

Direct comparison was made between experimental and numerical thmst 
measurements. As expected, the panel code predicts thmst consistently higher than that 
achieved in the experiments. The difference between numerical and experimental thmst 
may be attributed primarily to separation, three-dimensional effects, and wing flex in the 
experimental model. In general, the experimental thmst values follow the numerical 
value trends. 
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V. RECOMMENDATIONS 


The parameter space associated with flapping-wing propulsion is broad. An 
infinite variety of flapping-wing configurations and dynamics are available for analysis. 
The difficulty of evaluating this space is greatly eased by numerical methods. Programs 
can be easily written to continuously change and analyze potential candidates. A 
program for analyzing output data would contribute to the speed of analysis. The 
numerical analysis can play a vital role in determining promising experimental 
configurations. 

Experimentally, efforts should be greatly eased by the tunnel fan overhaul and 
general clean-up, installation of new equipment for data collection, standardization of 
measurement techniques, and standardization of data reduction. Promising 
configurations can be readily tested. Experimental testing could be improved with the 
addition of a sting to measure lift and thrust simultaneously. 

Newly acquired smoke equipment and high-speed photographic equipment should 
greatly aid in smoke visualization studies. Flow visualization will provide a better 
understanding of flow dynamics, primarily by identifying flow separation. In addition, 
the installation of the LDV system and LDV probe traverse mechanism should allow for 
mapping the complex flow field of the flapping wings. 

The identification of the source of divergence between numerical and 
experimental results could be investigated. The difference is currently attributed to 
separation, three-dimensional effects, and wing flex in the experimental model. 
Experimentally, 3-D effects can be minimized with tip plates, either attached to a trailing 
stationary wing or external to the model. Wing tip flex and phase angle can be measured 
for application to the numerical code. A numerical strip analysis can be done by applying 
the model’s varying spanwise wing deflection and phase angle to the code. 
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APPENDIX A. EXPERIMENTAL EQUIPMENT FAILURES 


The scope of this thesis has been limited by time spent repairing old equipment 
and installing new equipment. The initial experimental effort was delayed by a tunnel 
fan failure. Fan disassembly indicated worn bearings. The fan assembly was removed, 
cleaned, bead blasted, and the blades painted. The 32 fan bearings, one thrust and one 
roller bearing for each of the 16 blades, were removed and replaced. The tunnel motor 
was inspected and motor bearings greased. Fan shaft bearings were inspected and 
greased. The fan was reassembled and blades reset with uniform reduced pitch on each 
blade. Check runs following repairs and subsequent experimental work has demonstrated 
significantly reduced vibration levels. The reduction in blade pitch has changed the 
tunnel velocity range from 1 - 10.5 m/s to 0 - 9.5 m/s. 

Initial LDV velocity measurements following tunnel repair indicated high vertical 
velocity at the test section. It was suspected that tunnel fan swirl was effecting upstream 
flow. Adjustable vanes just upstream of the fan were in a state of disrepair. The vanes 
were removed and honeycomb installed. In addition, general tunnel clean-up included 
removal of old experimental fixtures, tunnel floor refmishing, and smoothing of tunnel 
walls and ceiling by sealing holes and smoothing seams. The old pitot-static tube and the 
associated fluid manometer were removed and replaced by a better positioned pitot-static 
tube with a differential pressure transducer. A flood damaged 3-axis traverse mechanism 
was disassembled and cleaned, one motor replaced, and mounted at the test section 
window to support LDV measurements. 

Subsequent to the installation of downstream honeycomb and tunnel clean-up, 
LDV velocity measurements again indicated high vertical velocities. An examination of 
the LDV data indicated an inconsistent output. A series of difficulties were experienced 
with the LDV system, culminating with loss of lasing. A technician realigned and 
cleaned the laser mirrors. The old LDV processing equipment was replaced by new 
equipment. The signal processor was replaced (IFA 755 for EFA 550), and a Colorlink 
Plus Receiver installed. A new DMA board and software were installed on a new 
computer. The new LDV system has posed no problems. A description of the LDV 
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system can be found in the Experimental Apparatus section and Appendix C. Tunnel 
turbulence intensity levels are discussed in Appendix D. 

The flapping-wing mechanism’s motor armature was burned up during one 
experimental run. The motor was removed and replaced. The new motor’s voltage was 
limited to 26V for a maximum of one minute. This voltage limit restricted maximum 
flapping frequency to 7 Hz. 

The smoke equipment proved difficult to use for LDV measurements. Recently 
acquired smoke equipment and high-speed photographic equipment should greatly aid in 
smoke use for LDV and visualization studies. 

The accuracy of the oscilloscope used for voltage measurements came into 
question during experimental runs. The oscilloscope was replaced with a Link 
Instmments DSO-2102 digital storage oscilloscope. The oscilloscope connects to a 
desktop PC for display and recording of the voltage signals. The oscilloscope is 
described in the Miscellaneous Equipment section. 

Experimental configurations can now be more readily tested. The tunnel fan 
overhaul and general clean-up, installation of new equipment for data collection, 
standardization of measurement techniques, and standardization of data reduction should 
greatly ease experimental efforts. 



Figure 20. Tunnel Fan Complete. 
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APPENDIX B. NUMERICAL DATA 


Two NACA 0014 airfoils, one atop the other, equal chordlength, cCgeom=0, separated by 
mean of 1.4c, plunge-only motion w/ plunge amplitude h/c=0.4 


USPOT Data 


USPOT at 2.95 Hz 


Vel 

Thrust 

Eff 

Red freq 

[m/s] 

[N] 

Vt 

k 

11.859 

0.08058 

0.92475 

0.10 

5.929 

0.07534 

0.93556 

0.20 

4.941 

0.07160 

0.92392 

0.24 

4.235 

0.06768 

0.90946 

0.28 

3.706 

0.06381 

0.89338 

0.32 

3.294 

0.06011 

0.87642 

0.36 

2.965 

0.05664 

0.85910 

0.40 

2.695 

0.05346 

0.84177 

0.44 

2.372 

0.04921 

0.81630 

0.50 

1.976 

0.04333 

0.77646 

0.60 

1.482 

0.03567 

0.70877 

0.80 

1.186 

0.03128 

0.65492 

1.00 

0.847 

0.02709 

0.57124 

1.40 

0.593 

0.02503 

0.47743 

2.00 


USPOT at 
Vel 

7.03 Hz 
Thrust 

Eff 

Red freq 

[m/s] 

[N] 

Vt 

k 

28.273 

0.45804 

0.92475 

0.10 

14.137 

0.42827 

0.93556 

0.20 

11.781 

0.40699 

0.92392 

0.24 

10.098 

0.38475 

0.90946 

0.28 

8.835 

0.36273 

0.89338 

0.32 

7.854 

0.34167 

0.87642 

0.36 

7.068 

0.32198 

0.85910 

0.40 

6.426 

0.30387 

0.84177 

0.44 

5.655 

0.27975 

0.81630 

0.50 

4.712 

0.24633 

0.77646 

0.60 

3.534 

0.20276 

0.70877 

0.80 

2.827 

0.17778 

0.65492 

1.00 

2.020 

0.15396 

0.57124 

1.40 

1.414 

0.14226 

0.47743 

2.00 


USPOT at 5.01 Hz 


Vel 

Thrust 

Eff 

Red freq 

[m/s] 

[N] 

Vt 

k 

20.138 

0.23238 

0.92475 

0.10 

10.069 

0.21727 

0.93556 

0.20 

8.391 

0.20648 

0.92392 

0.24 

7.192 

0.19520 

0.90946 

0.28 

6.293 

0.18403 

0.89338 

0.32 

5.594 

0.17334 

0.87642 

0.36 

5.035 

0.16335 

0.85910 

0.40 

4.577 

0.15416 

0.84177 

0.44 

4.028 

0.14193 

0.81630 

0.50 

3.356 

0.12497 

0.77646 

0.60 

2.517 

0.10287 

0.70877 

0.80 

2.014 

0.09019 

0.65492 

1.00 

1.438 

0.07811 

0.57124 

1.40 

1.007 

0.07217 

0.47743 

2.00 
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APPENDIX C. LASER DOPPLER VELOCIMETRY (LDV) 

A. INTRODUCTION 

Laser Doppler Velocimetry (LDV), also referred to as Laser Velocimetry (LV), 
was utilized for measuring flow velocity and fluctuations. Although sophisticated and 
expensive, it provides a non-intrusive, very accurate means of measuring fluid particle 
velocity. Its dependence on optical access and particles of proper size in the flow to 
measure velocity can introduce complications. A general description is provided under 
Experimental Apparatus, the details of LDV position and sampling under Experimental 
Data Collection, and LDV use for tunnel calibration under a separate Appendix. This 
section is intended to provide more details of the system used and the concepts associated 
with LDV. 

B. LDV APPARATUS 

The LDV system is depicted in Figure 22. The system uses a 5W, water-cooled 
Coherent Innova Series Argon-ion laser for the light source. A TSI Model 9201 
Colorburst multicolor beam separator unit is used to separate the laser colors (green - 
514.5nm, blue - 488nm, and purple - 476nm), split the beams and frequency shift one 
beam in each color. After splitting the beams, one of the two blue beams and one of the 
two green beams was frequency shifted 40 MHz. The purple channel was blocked. The 
beams were coupled into a 2-component, 4-beam, Model 9832 fiber-optic probe through 
a Model 9271 coupler (one for each beam). The beams were transmitted to the probe via 
optical receiving fiber. Probe beam spacing was 50 mm and the lens focal length was 
350 mm. The photo of the probe at the test section window is shown in Figure 23. A 
multi-mode receiving fiber located in the probe received and forwarded the reflected light 
to a Model 9160 photo multiplier tube (PMT). Although the probe is a two-component 
unit, only one component of the receiving optics was available. Thus, only single 
component measurements were possible at a given time. The output of the PMT was 
connected to a 9186A frequency downshifter. The downshifter is connected to a TSI 
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Model IFA-755 signal processor via the Colorlink Plus Multicolor Receiver. The 
Colorlink Plus contains an integral PMT and frequency downshifter. The probe 
connector in use, however, is incompatible with the Colorlink Plus. As such, an older 
stand-alone PMT and downshifter were used. In this setup, the Colorlink Plus serves no 
purpose except to relay the signal from the downshifter to the rFA-755 signal processor. 
The signal processor distinguishes the signal burst from the noise based upon the signal- 
to-noise ratio. The processor determines frequency, phase, burst transit time, and burst 
arrival time and sends these signals to the computer. The DMA board and FIND for 
Windows software package on the computer converts the signals into velocity, and is 
used for data display and storage. 

LDV requires no calibration, but does require a simple alignment. The probe was 
aimed at a black wall. Each of the three colors in the Color Burst Module have 
adjustments for changing the beam’s X, Y and Z position as well as tilt. Adjustments 
were made to the blue and the green laser spots to obtain the highest intensity (brightest). 
These adjustments ensure that the beam entering the fiber-optic cable is within the 
Brewster angle limitations. Significant signal losses occur if the beam does not enter the 
fiber-optic cable at the proper angle. 

C. LDV PROCESSING CONCEPTS 

1. Interference fringe spacing. 

Interference fringes are formed when two light waves interfere at a crossing point. 
Interference fringes are a set of bright and dark bands across the probe volume, formed in 
the plane of the beams. The component of velocity measured is perpendicular to the 
bisector of the angle formed by the beams, i.e. perpendicular to the fringes. The fringe 
spacing value is required for velocity calculations. The FIND software determines fringe 
spacing automatically with either the focal length and beam spacing, or the half angle, K. 
A pictorial representation of the LDV beam geometry is found in Figure 21 and a 
summary of the calculations involved follows. 
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Fringe spacing, df, is related to the wavelength of the incident light, X, and the 
half angle between the two beams, K, by: 


LL f — - 

2sinK' 

The frequency detected at photodetector,can be 
the x-direction, Vx by the following: 


(14) 

related to the free-stream velocity in 


fD=^ (15) 

df 

Trigonometry establishes the relationship between the beam spacing, d, and the focal 
length,/. 


d 

sin K - , = 


2 / 

cosy =-7 

44^77 


(16, 17) 


For the probe used, the beam spacing is fixed at 50 mm, and a lens with a focal length of 
350 mm was used exclusively. 

Combining Equations 14 and 16 gives: 



A-^4f-+d^ 

2d 


(18) 


This Equation demonstrates that fringe spacing is constant for a given LDV 
configuration. Calculating the fringe spacing from the known values of the wavelength 
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of green light, focal length, and beam spacing yields a fringe spacing of 3.6107 |xm. The 
half angle equals 4.086°. [Ref. 15] 


2. Signal Processor 

The IFA 755 signal processor has a low-pass and a high-pass filter to limit data 
collection to Doppler frequencies of interest. The tunnel velocity range is approximately 
0 tolO m/s. Applying Equation 15, the expected Doppler shift frequency,^, was 
calculated. The low and high-pass filter settings for the velocities of interest are 
contained in Table 2. 


LDV Low and Hi-Pass Filter 

Veloci 

ty Range 

Frequenc; 

y Range 

Green 

Blue 





df=:3.61C 

17 mm 

df=3.4247 mm 

Low 

Low [Hz] 

Hi 

Hi [Hz] 

Low [m/s] 

Hi [m/s] 

Low [m/s] 

Hi [m/s] 

30 kHz 

30000 

300 kHz 

300000 

0.108 


■kies 

1.027 

100 kHz 

100000 

IMHz 

1000000 

0.361 


■BBS 

3.425 

300 kHz 

300000 

3MHz 

3000000 

1.083 

10.832 

1.027 

10.274 


Table 2. LDV Doppler Frequency Filter Settings. 

The colorburst multicolor beam separator unit was used to split the blue and green 
beams and frequency shift one beam of each color. The split beams were frequency 
shifted 40 MHz. The frequency downshifter shifted this frequency 39.8 MHz, for an 
effective frequency shift of 0.2 MHz (200 kHz). As the frequency goes down, the system 
resolution goes up. Downshifting, therefore, improves resolution to allow more precise 
measurement of the flow velocity. A zero velocity flow then would have a frequency of 
200 kHz measured at the photo detector (i.e. no Doppler shift). Positive flows should 
have frequencies above 200 kHz and negative flows less than 200 kHz. 


3. 

The 


Probe Volume Characteristics 

properties of the probe volume are calculated by the following equations: 

A2f dd 

d.3 =-^, (19,20,21) 

* nD.. cos a: sinx" 
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(22, 23) 


TTd^.. 

Vol = -- 

6cos/£:sinK' 



l.27d 

Z)., 

e " 


where dg is the diameter of the intersecting region of the probe volume, Dg is the diameter 
of the beam at the focusing lens (De=l.O mm), d^ is the probe volume diameter across the 
fringes, is the probe volume length, Vol is the probe volume, and Nfr is the number of 
fringes. A summary of probe volume characteristics, including the application of the 


above equations, is listed in Table 3. 


Characteristic 

Value 

Color 

Green 

Blue 

Wavelength, X 

514.5 nm 

488.0 nm 

Lens focal length,/ 

350 mm 

Beam spacing, d 

50 mm 

Half angle between beams, k 

4.086“ 

Fringe Spacing, df 

3.6107 |j,m 

3.4247 ^im 

Diameter of beam at focusing lens, Dg 

1.0 mm 

Beam diameter, d 

e 

229.28 nm 

217.47 urn 

Probe volume diameter, dm 

229.86 |im 

218.02 jxm 

Probe volume length, Im 

3.218 mm 

3.052 mm 

Probe volume, Vol 

0.0888 mm" 

0.0758 mm^ 

Number of Fringes, Nrr 

63 

.7 


Table 3. Summary of Probe Volume Characteristics. 


4. Velocity and Turbulence Intensity 

The Flow Information Display (FIND) software calculates the mean velocity of 
the flow by taking the mean of the samples taken at one data point. The number of 
samples taken per data point is selected by the user. The software calculates the velocity 
by multiplying the Doppler sift by the fringe spacing. That is. 
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Fo=f,+U=f,+^ 


(24) 


whrere Fd is the Doppler frequency measured at the photo detector,,^ is the selected 
frequency shift (200 kHz downshift),^ is the Doppler shift caused by flow velocity, and 
V is the speed of the flow (0 - 9.5 m/s). The software also determines the turbulence 
intensity at each data point. 

5. Seeding and Data Rate 

LDV is an optical method measuring the frequency shift of laser light scattered by 
particles passing through the probe volume. The probe volume is an ellipsoid formed by 
the crossing of the laser beams. LDV measures the velocity of the particles passing 
through the probe volume, not the fluid. If no particle passes through the volume, no 
velocity is measured. System accuracy is not only dependent upon the presence of 
particles, but also upon the particles following the flow, and on the particles scattering 
light well. Particles that meet these criteria are approximately 1 pm. If the particles are 
too small, they scatter light that is not strong enough for a good signal but are strong 
enough to be accepted by the processors. The result is a noisy signal. In most 
experiments with air, artificial seeding is required. Flow seeding is done by atomizers 
that produce uniform particles that disperse in the flow. A Rosco Model 4500 smoke 
generator was used in this experiment. The data gathered produced a statistically well- 
defined scatter of velocity. Improper seeding (either undersize particles or overseeding) 
would have resulted in noise. 


50 



DMA Cable 



Figure 22. LDV Setup. After Ref.s 17 and 18. 
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APPENDIX D. WIND TUNNEL CALIBRATION 


Following repair and refurbishment of the wind tunnel, a tunnel calibration was 
conducted using the LDV system and pitot-static tube. A general description of the LDV 
and pitot-static system can be found in the Experimental Apparatus section. A detailed 
description of the LDV is found in Appendix C. 

In preparation for calibration, the model and other associated equipment were 
removed from the tunnel. Measurements were taken over the tunnel velocity range of 0 
to 9.5 m/s by both the pitot-static tube and the LDV probe. The data was used to plot a 
curve comparing the LDV to pitot measurements. Turbulence intensity was recorded by 
the LDV for plotting as well. The last recorded evaluation of the turbulence intensity was 
performed by Costello [Ref. 19] with hot-wire anemometers. In the velocity range of 1.8 
to 9.8 m/s, he recorded a turbulence intensity range from 0.47% to 0.97%. A 750 mm 
focal length LDV lens will be acquired soon, allowing LDV measurements to the center 
of the tunnel, and an expanded turbulence intensity evaluation. 

The newly installed pitot tube is located approximately 2.5 m upstream of the 
model and LDV probe. The tube is 0.75 m above the tunnel floor, and can be traversed 
in and out of the flow, from 0 to 0.5 m across the tunnel. Further downstream, the LDV 
is mounted outside the tunnel on a 3-axis traverse mechanism abeam the model’s usual 
position, labeled as test section in Figure 13. 

An MKS Baratron type 223B differential pressure transducer connected to the 
pitot-static probe. The transducer voltage output is linear with pressure. Bernoulli’s 
equation was used to relate the pitot-static pressure to velocity: 


V = 


(25) 


where Vp,,^ is flow velocity as measured by the pitot-static tube, AP is the pitot-static 
pressure differential, and p is atmospheric density. 

For the first experiment, the LDV was moved from its normal position at the test 
section window to a plexiglass window that the pitot tube is mounted in. The two were 
essentially co-located, with the LDV probe 0.1 m above and 0.05 m in front of the pitot 
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tube. Both were positioned to take velocity measurements approximately 0.35 m from 
the tunnel wall, the range limit of the current LDV probe lens. LDV and pitot-static 
measurements were simultaneously taken and recorded for comparison. Pitot voltages 
were recorded by a digital voltmeter for 60 seconds, the time found to allow for 
consistent velocity convergence. The LDV measurement runs consisted of 5,000 samples 
taken over a period of 5 to 15 seconds, dependent upon smoke density. Pitot voltages and 
LDV free-stream velocity and associated velocity fluctuations were recorded for 
subsequent data reduction. 

Three curves were produced. The LDV velocity as a function of pitot-static 
voltage output is shown in Figure 24. Figure 25 is a linear fit of LDV and pitot velocity 
using the experimental data fit equation from Figure 24. Figure 26 plots the tunnel 
turbulence intensity versus tunnel velocity as recorded by the LDV. The turbulence level 
remained below 3.3% above 2 m/s. Tabular data is attached at the end of the Appendix. 
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Figure 24. LDV Velocity vs Pitot Voltage (LDV, Pitot Co-located). 
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Pitot Velocity [m/s] 


















In the second experiment, the LDV probe was mounted on its traverse, 2.4 m 
downstream from the pitot-static tube. The pitot-static tube remained at its fixed height 
of 0.75 m, with the LDV 0.1 m lower. Both were positioned to obtain the velocity 0.35 
m from the tunnel wall. The tunnel cross-section is nominally 1.5 x 1.5 m, but does 
expand 2.6% over the 2.4 m between the pitot-static tube and LDV probe. It is assumed 
this was designed in to account for boundary layer growth. 

The first experiment’s procedures were repeated with the LDV at the test section. 
LDV velocities were lower than the pitot-static velocities, as would be expected due to 
the expanding area. Figure 27 shows the LDV velocity as a function of pitot-static 
voltage. Figure 28 is a near-linear fit of LDV and pitot velocity using the equation from 
Figure 27. Figure 29 presents the LDV recorded turbulence intensity at the position 
normally occupied by the model. The turbulence level remained below 1.75% above 1.5 
m/s. In general, the intensity levels at the two tunnel positions were comparable. Tabular 
data follows the figures. 



Pitot Voltage [mV] 


Figure 27. LDV Velocity vs Pitot Voltage (LDV downstream of pitot at Test Section). 
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Turbulence Intensity [%] 


Pitot Velocity [m/s] 



LDV Mean Velocity [m/s] 

Figure 29. Turbulence Intensity vs LDV Velocity at Test Section 


















Pitot/LDV Co-located 








Run Sample LDV Mean 

Vel 

Turbulence 

Turbulence 

Pitot 

Pitot 

Pitot Vel 

Pitot Vel 

Size Velocity 

Std Dev 

Intensity 

Intensity 

(60 s mean) 

(60 s mean) 

fm Eqn 1 

fm Eqn 2 

[m/s] 

[m/s] 

[%] 

[m/s] 

[mV],[V] 

[mV] 

[m/s] 

[m/s] 


#NUM! 

0.534" 


#NUM! 

0.552 


.029 

.510 


2.012 

2.489 






18 _ 5000 8.3 

19 _ 5000 8 . 8 < 

20 5000 9.2! 

21 I 50001 9.4: 

Eqn 1=14.7r(volt$)^5 
Eqn 2=0.4819{mV)^.489 


93.437 

115.434" 


140.493 
165.057 " 
195.931 " 
0.227 " 
0.258 " 
0.295 " 
0.336 ~ 
0.374" 
0.4171 
0.429’ 


93.437 

115.434 


140.493 

165.057 

195.931 

226.8981 

258.016 

295.475 

336.066 

374.441 

416.852 

428.836 


Eqn 1 based on Bernoulli | 

Eqn 2 derived fm experimental data 


Table 4. Velocity Data for Pitot-LDV Co-located. 


LDV on Traverse (downstream 
Run I Sample iLDVMeanl 


Turbulence | 


Intensity 

[%] 




Turbulence 


Intensity 

[nn/s] 


Pitot 


(60 s mean) 
[mV], [V] 


Pitot 


(60 s mean) 
[mV] 


>0.14 

1.36' 


Pitot Vel I 


fm Eqn 1 
[m/s] 


#NUM! 

0.543 


Pitot Vel 


fm Eqn 2 
[m/s] 


#NUM! 

0.518 



Eqn 1=14.71 *(volts)^5 
Eqn 2=0.443(mV)^.5045 


_ 1.330 _ OA 

_ 1.457 _ OS 

1.167 0.( 


.218 _ OA 

.725 0.( 

I 0.901 _ OA 

.290 _ OX 

.224 _ OX 

1.968 _0^ 

.866 _ 0 ^ 

.166 _ OX 

.624 _ OX 

.941 0.( 

.353 0.- 

.958 _OX 

.057 _ OX 

I Eqn 1 based on Bernoulli _ 

[Eqn 2 derived fm experimental data 


29.49 
41.711 
56.87"" 


73.64 

93.10 ” 

115.46 

140.10| 

165.50 

195.04 

0.23 

0.26 

0.30 

0.34 

0.37 

0,42 

0.43 



Table 5. Velocity Data for Pitot-LDV with LDV at Test Section. 
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APPENDIX E. MODEL THRUSTAVEIGHT CALIBRATION 


Weiqht vs Voltaae Calibration 




















Weight [g],| 

[N] 

[g] 

0 

3 

5 

10 

20 

35 

50 

65 

75 

[N] 

0.0000 

0.0294 

0.0490 

0.0981 

0.1961 

0.3432 

0.4903 

0.6374 

0.7355 

Run a mean [V] 

0.0041 

0.0877 

0.1402 

0.2603 

0.5192 

0.8966 

1.2710 

1.6630 

1.9240 

Run a std dev [V] 

0.0014 

0.0005 

0.0006 

0.0010 

0.0007 

0.0006 

0.0006 

0.0007 

0.0008 

Run b mean [V] 

0.0047 

0.0851 

0.1383 

0.2539 

0.5168 

0,9026 

1.2901 

1.6682 

1.9351 

Run b std dev [V] 

0.0059 

0.0007 

0.0006 

0.0008 

0.0010 

0.0011 

0.0015 

0.0010 

0.0011 

Run c mean [V] 

0.0060 

0.0875 

0.1389 

0.2575 

0-5169 

0.8951 

1.2862 

1.6686 

1.9203 

Run c std dev [V] 

0.0043 

0.0012 

0.0009 

0.0009 

0.0009 

0.0010 

0.0008 

0.0007 

0.0006 

Avg Voltage [V] 

0.0049 

0.0868 

0.1391 

0.2572 

0.5176 

0.8981 

1.2824 

1.6666 

1.9265 

Avg Distance [mm] 

DEESES 


1.3914 

2.5723 

5.1764 

8.9810 

12.8244 

16.6659 

19.2649 

Avg Std Dev [V] 



0.0007 


0.0009 

0.0009 

0.0010 

0.0008 

0.0008 

Corrected Voltage [V] 

0.0000 


0.1342 

0.2523 

0.5127 

0.8932 

1.2775 

1.6616 

1.9215 

Thrust fm eqn 

0.0000 

0.03141 0.0514 

0.0967 

0.1965 

0.3423 

0.4897 

0.6369 

0.7365 

Error [%] 

#DIV/0! 

6.61351 4.9034 

1.3925 

0.1951 

0.2578 

0.1366 

0.0820 

0.1398 
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APPENDIX F. MATLAB PROGRAM FOR DSO-2102 OSCILLOSCOPE DATA 

REDUCTION 


% LT T.C. Lund 

% DSO-2102 PC Oscilloscope data conversion program 
% Filename: Convert.m 


% This program converts data saved from the DSO-2102 oscilloscope 
% software into voltage, range, thrust, and frequency. It converts 
% and saves only the first two columns of the .prn file type, but 
% could easily be modified to read all six columns, or to read the 
% .csb or .CSV data file types. 


% For program users: 

% - When you take your data, save your .prn files as a filename of 
% less than eight characters total, with the last characters as 
% sequential numbers (e.g. thrstlOO.prn, thrstlOl.prn, etc). 

% “ Enter the data under "Fill in this part." 

% - Output goes to filename.wri (e.g. thrstlOO.wri, thrstl01.wri) and 
% Mean_Thrust.wri (Caution: The program deletes existing 
% filename.wri and Mean_Thrust.wri files). The programs open in 
% Wordpad for quick viewing, but the cols will be out of alignment. 

% The output is designed for Excel, where the cols will be aligned. 

% - You also must update the following eqns with each configuration: 

% 1. The 5 eqns that convert voltage to thrust. 

% 2. The 1 eqn for weight calibration error. 

% - Program assumes you select +/-5V on laser rangefinder (not O-IOV) 

% - Program assumes range data goes into DSO Chan B, and trip switch 
% goes into DSO Chan A. 

% - Program assumes that the frequency trip switch is a negative 
% voltage. See notes under DSO data variables and 3rd if statement, 
clear 

delete Mean_Thrust.wri; 


filename='thrst'; 
filetype='.prn'; 
startfilenumber=100; 
lastfilenumber=163; 
voltsperdiv=.2; 
sample_rate=1000; 
v_zero=.00251; 


Fill in this part *★*****★**★★★***★**★*★★*★***** 


% First file number 
% Last file number 
% Volts per division fm DSO 
% Sample rate in Hz fm DSO 
% DSO zero voltage offset 


fidl = fopen('Mean_Thrust.wri','w'); 

fprintf(fidl,'Filename\tVel\tVel\tMean\tCorrectedXtTotal Error\t'); . 
fprintf(fidl, 'Avg Dev\tLaser\tDSO\tWt Cal\tMean\tAvg Dev\tMin\tMax\n'); 
fprintf(fidl,'\t\tStd Dev\tThrust\tThrust\t+/-XtErrorXtErrorVtErrorNt'); 
fprintf (fidl, 'Error\tFrequency\tFrequencyXtFrequencyXtFrequency\n') ; 
fprintf (fidl, '\t[m/s] \t [m/s] \t [N] \t [N] \t[N] \t [N] \t[N] \t [N] \t[N] \ t'); 
fprintf(fidl, ' [Hz]\t[Hz]\t[Hz]\t[Hz]\n\n'); 
fclose(fidl); 
for i=100:lastfilenumber 

clear data % Clears data matrix each iteration 
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% ********************** Loads data ***’^******★★★*★****★★******★★*** 

disp( ['* 1 ^***** Computing ' , filename,num2str(i) , filetype, ' **★*★**']) 

rawdata=load ([filename,num2str(i),filetype]); 

data(:,1)=rawdata{:,2); % CH B is loaded to data col 1 

data{:,8)=rawdata{:,1); % CH A is loaded to data col 8 

% **★★★★★* Converts data to voltage and distance (lV=lcm) ********** 

data(;,2)={voltsperdiv*5.12)-((data(:,!))*(voltsperdiv^lO.24/256)); 

data(:,9)=(voltsperdiv*5.12)-((data(:,8))*(voltsperdiv^lO.24/256)); 

% ★**★★*** Voltage correction ************************************** 
data{:,3)=data(:,2)+v_zero; % DSO volts corrected for Vzero offset 
% ★★*****★ Converts voltage to thrust [N] ****★★★**★*★★★★★★★★****★** 
data(:,5)= 0.3833* data(:,3); 

% *★***★*★ Computes the avg dev of thrust & freq in data *★★★*★***** 

% DSO data variables 

sig_trip = 175; % The minimum data value that trips micro-switch 

sig_oldl = 257; % sig_oldl, 2, and 3 are the 3 data points prior 

sig_old2 = 257; % to the current pt. Use 257 if the trip signal 

sig_old3 = 257; % is neg voltage. Use -1 if trip is a + voltage. 

% Voltage variables 

cycle_volts=0; % Sum of voltage for one cycle 

cycle_mean_volts=0; % Mean volts of one cycle 
volts_siim=0; % Sum of voltages of all cycles 

% Thrust variables 

cycle_thrust=0; % Sum of thrust for one cycle 

cycle_mean_thrust=0; % Mean thrust of one cycle 

thrust_sum=0; % Sum of mean thrust of all cycles 

% Frequency variables 

cycle_mean_freq=0; .% Mean frequency of one cycle 

freq_sum=0; % Sum of frequencies of all cycles 

freq=0; % Variable for min frequency 

k=l; % Min freq counter 

% Counter variables 

nr_cycles=“l; % Counter for total number of cycles in data 

cycle_data_pts = 0; % Number of data pts within one cycle 

total_data_pts=0; % Number of data pts of all cycles 

data^oints = length(data ( :, 1) ) ; % Counts rows of data matrix 
for j = l :data_points % Looks at all data points 
if abs(data{j,2)) > 5 

error ('** Voltage is outside rangefinder limits {+/- 5V) **') 

end 

if nr_cycles >= 0 % Ignores partial data before 1st cycle 

cycle_volts = cycle_volts + data (j-1,2); % Sums cycle volts 
cycle__thrust = cycle_thrust + data(j-l,5); % Sums cycle thrust 
cycle_data_pts=cycle_data_pts+l; % Counts data pts in cycle 

end 

% Is micro-switch tripped? 

if data(j,8)>=sig_trip & sig_oldl<sig_trip & ... 

sig_old2<sig_trip & sig_old3<sig_trip % This assumes a neg 
% trip voltage. Change >, < signs if trip is a + voltage. 
nr_cycles=nr_cycles+l; % Count only when ifs met 

if nr_cycles > 0 % Ignores partial data before 1st cycle 
% Volts 

cycle_mean_volts=cycle_volts/cycle_data_pts; % Cycle mean 
data(j-1,4)=cycle_mean_volts; % Stores cycle mean volts 
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volts_siiin=volts_siiin+cycle_mean_volts; % Sums total volts 
% Thrust 

cycle_mean_thrust=cycle_thrust/cycle_data_pts; % Cycle mean 
data(j-l,6)=cycle_mean_thrust; % Stores cycle mean thrust 
thrust_sum=thrust_sum+cycle_mean_thrust; % Sums tot thrust 
% Frequency 

cycle_mean_freq=sample__rate/cycle_data_pts;. % Cycle mean 
data(j-l,7)=cycle_mean__freq; % Stores cycle mean frequency 
freq_sum=freq_sum+cycle_mean_freq; % Sums all cycle freqs 
if data(j“1,7)>0 

freq(k) = cycle_mean_freq; % Save nonzero freqs 
k=k+l; 

end 

total_data_pts=total_data_pts+cycle_data_j5ts; 

% Reset variables 

cycle_volts=0; % Reset cycle voltage for new cycle 

cycle_thrust=0; % Reset thrust summation for new cycle 
cycle_freq=0; % Reset freq summation for new cycle 

cycle_data_pts=0; % Reset nr of data pts for new cycle 

end 

end 

sig_oldl = data(j,8); % Last signal value 

if j>=3 

sig_old2 = data(j-l,8); % The data value before the last one 
sig_old3 = data(j-2,8); % Three data points back 

end 

end 

% Voltage and frequency calculations: 
volts_mean=volts_sum/nr_cycles; 

if volts_mean==0 % No or very slow flap case, still need data 

volts_mean=mean(data(:,3)); 

end 

inner_suml=0; % Inner summation for adev of volts 
inner_sum2=0; % Inner summation for adev of freq 

freq_mean=freq_sum/nr_cycles; % Mean freq of cycles 
max_freq=max(freq); % Max cycle freq of all cycles 
min_freq=min(freq); % Min cycle freq of all cycles 
for m=l:data_points 
if data(m,6)>0 

inner_suml=inner_suml + abs(data(m,4)-volts_mean); 
inner_sum2=inner_sum2 + abs(data(m,7)-freq_mean); 

end 

end 

adev_freq=inner_s\am2/nr_cycles; % Avg dev of frequency 
% Thrust 

thrust_mean=0.3833*volts_mean; % Mean thrust of cycles 
% *** There are 4 sources of thrust error: 

% *** 1st one is due to cyclic nature of flapping-wing thrust: ***** 
adev_volts=inner_suml/nr_cycles; % adev fm voltage data 
% *** 2nd source of error - laser rangefinder error [mm] *********** 
if abs(volts_mean) < 3.5 

error_laser = 0.1 + 0.002 * abs(10 * volts_mean); 
else 

error_laser = 0.25 + 0.005 * abs(10 * volts_mean); 
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end 

error_laser = error_laser/10; % Convert range error to volts error 
% *** Oscilloscope volt max error=resolution of DSO-2102-3rd error * 
error_DSO = (voltsperdiv’^S . 12)/128 ; 

% *** Convert errors fm voltage to thrust [N] error **************** 
adev_thrustl=0.3833*adev_volts; 
error__thrust2 = 0.3833*error_laser; 
error_thrust3=0.3833*error_DSO; 

% *** Computes adev due to weight calibration - 4th error: ********* 
thrust_hi= 0.3836*volts_mean - 3E-06; 
error_thrust4 = thrust_hi-thrust__mean; 

error_thrust=adev_thrustl+error_thrust2+error_thrust3+error_thrust4 ; 

% **★★★★**★******★ Saves all data to filename.wri ****************** 

filesave2=([filename,num2str(i),'.wri']); 

fid2 = fopen(filesave2,'w'); 

fprintf{fid2,[filename,num2str(i)]) ; 

fprintf(fid2,'\nCH B\tCH B\tCH B\tCH B\tCH B\t'); 

fprintf(fid2,'CH B\tCH B\tCH A\tCH A\n'); 

fprintf(fid2,'Data\tVoltage\tAbsol VoltXtCycle Mean\tThrust\t'); 
fprintf(fid2,'Cycle MeanXtFreq MeanXtDataXtVoltageXn'); 
fprintf(fid2,'Xt[V], [cm]Xt[V], [cm]Xt[V]Xt[N]Xt[N] Xt'); 
fprintf(fid2,'[Hz]XtXt[V]Xn'); 

fprintf(fid2,'XnXtMean thrust =XtXtXt %10.8f',thrust_mean); 
fprintf(fid2,'XnXtThrust Avg Dev Error =XtXtXt %10.8f',adev_thrustl); 
fprintf (fid2 ,'XnXtLaser Error =:XtXtXt %10.8f ' , error__thrust2) ; 
fprintf{fid2,'XnXtOscilloscope Error =XtXtXt %10.8f',error__thrust3); 
fprintf(fid2,'XnXtWt/Thrust Cal Error =XtXtXt %10.8f',error_thrust4); 
fprintf(fid2,'XnXtTotal Error =XtXtXt %10.8f',error_thrust); 
fprintf(fid2,'XnXtXtXtXtMean frequency =XtXt %4.2fXn',... 
fr’eq^.mean) ; 

fprintf(fid2XtXtXtXtAvg Dev =XtXt %4.2fXn',adev_freq); 
fprintf(fid2XtXtXtXtMin Freq =XtXt %4.2fXn',min_freq); 
fprintf(fid2,'XtXtXtXtMax Freq =XtXt %4.2fXnXn',max_freq); 
fprintf(fid2,'%3.OfXt %4.3fXt %5.4fXt %5.4fXt %8.5fXt %8.5fXt %5.2fXt 
%3.0fXt %4.3fXn',data.'); 
fclose(fid2); 

% ***★★**★★***★ Saves summary of data to Mean_Thrust.wri *********** 
fid3 = fopen('Mean_Thrust.wri','a'); 
fprintf(fid3,[filename,num2str(i)]); 

fprintf{fid3,'XtXtXt%10.8fXtXt%10.8fXt%10.8fXt%10.8fXt%10.8fXt 
%10.8fXt%4.2fXt%4.2fXt%4.2fXt%4.2fXn',... 

thrust_mean,error_thrust,adev_thrustl, error_thrust2,error_thrust3, . . . 
error_thrust4, frequmean, adev_freq,min_f req,max__freq) ; 
fclose(fid3); 

end 

disp(['********* Program Complete ***********']) 
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APPENDIX G. EXPERIMENTAL DATA 
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